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Abstract 

An improved version of the Olami-Feder-Christensen model has been intro- 
duced to consider avalanche size differences. Our model well demonstrates the 
power-law behavior and finite size scaling of avalanche size distribution in any 
range of the adding parameter p a dd of the model. The probability density functions 
(PDFs) for the avalanche size differences at consecutive time steps (defined as re- 
turns) appear to be well approached, in the thermodynamic limit, by g-Gaussian 
shape with appropriate q values which can be obtained a priori from the avalanche 
size exponent r. For the small system sizes, however, return distributions are found 
to be consistent with the crossover formulas proposed recently in Tsallis and Tir- 
nakli, J. Phys.: Conf. Ser. 201, 012001 (2010). Our results strengthen recent 
findings of Caruso et al. [Phys. Rev. E 75, 055101(R) (2007)] on the real earth- 
quake data which support the hypothesis that knowing the magnitude of previ- 
ous earthquakes does not make the magnitude of the next earthquake predictable. 
Moreover, the scaling relation of the waiting time distribution of the model has 
also been found. 
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1 Introduction 



Self-organized criticality (SOC) is a concept designed to describe extended dynamical 
systems reaching a statistically stationary state, characterized by power-law distribution 
functions in both space and time, without any "fine tuning" of an external parameter. 
SOC was first introduced as a subject by Bak, Tang and Wiesenfeld (BTW) in 1987 
JT]. In their well-known paper, they proposed a sandpile model and found the system 
showed SOC phenomenon with bulk conservation law and open boundary conditions 
H1EI- SOC has been proposed as a way to model the widespread occurrence of power 
laws, i.e., the abundance of long-range correlations in space and time in various sys- 
tems, such as chemical reactions, evolution, avalanches, forest burns, heart attacks, 
market crushes, earthquakes, etc J3l HJ . In order to forecast earthquakes, several sta- 
tistical models of earthquakes embodying such SOC features have been proposed and 
studied 16] |8] |9] 13 ■ F° r example, one is the Burridge-Knopoff (BK) model 0, in 
which an earthquake fault is modeled as an assembly of blocks mutually connected via 
elastic springs which are slowly driven by external force. Another extensively stud- 
ied statistical model might be the Olami-Feder-Christensen (OFC) model, which was 
first introduced by Olami, Feder and Christensen in 1992 as a simplification of the 
BK model. Mapping the BK model into a two-dimensional lattice, they simulated the 
earthquake behavior and introduced dissipation into the family of the SOC systems 
iflOl 151 [TP . Numerical studies have revealed that the OFC model exhibits apparently 
critical properties such as the Gutenberg-Richter(GR) law and the Omori law lfl2l . For 
these reasons, the OFC model has been regarded as a typical nonconservative model 
exhibiting SOC. 

Many works of OFC model have focused on the homogeneous lattice network 
lfT4l[T31[T3l . however, the actual transmission of seismic energy or force is often inho- 
mogeneous |[T6l[T7l . We know that earthquakes occur as a result of the relative motion 
of tectonic plates and the seismic energy will be released in the form of earthquake 
waves (primary wave or secondary wave). This process takes place from the epicenter, 
which is below the earth surface and spread through the elastic vibration of the rocks. 
Due to different geological conditions, the earthquake wave in the rock will spread with 
different velocities and rates of decay. This will cause different energy decay in differ- 
ent geological conditions, therefore the heterogeneity of energy transfer occurs. So it is 
reasonable to assume that the real earthquake system is heterogeneous, and people can 
easily conclude that the heterogeneous factor should be investigated in the earthquake 
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model. Recently, some works have already been carried out along these lines: Baiesi 
and Paczuski proposed a metric to quantify correlations between earthquakes based on 
scale-free networks. According to this metric, typical events are strongly correlated to 
only one or a few preceding ones ifTSl . Thus a classification of events as foreshocks, 
main shocks, or aftershocks emerges automatically. Epicenter network of OFC model 
has been investigated by Peixoto and Prado 11191 , in which they obtain a direct network 
and show a sharp difference between the conservative and nonconservative regimes. In 
the scale free and directional network models, the energy is released either randomly 
or uniformly. In contrast to them, our energy release relates to the nature of adjacent 
rocks. We notice that the tectonic plates which have higher stress are prone to be af- 
fected by other plates. It can collect more energy or force released by other plates. In 
order to simulate this phenomenon, we introduce edge weight which determines how 
the energy is transferred from one point to another in the coupled-map lattice, to inves- 
tigate the SOC behavior on the inhomogeneous network. This work aims to study the 
self-organized criticality behavior of the non-conservative improved OFC model. 



Original OFC model. In the OFC model, "stress" variable Fi (Fi > 0) is assigned to 
each site on a square lattice with L x L sites. Initially, a random value in the interval 
[0, 1] is assigned to each Fi, where Fi is increased with a constant rate uniformly over 
the lattice until, at a certain site i, the Fi value reaches the threshold, F t h = 1. Then, 
the site i "topples" and a fraction of stress aFi (0 < a < 0.25) is transmitted to each 
of its four nearest neighbors, while Fi itself is reset to zero, namely, 



where "nn" denotes the set of nearest-neighbor sites of i. If the stress of one "nn" site 
j exceeds the threshold, i.e., Fj enn > F t h = 1, the site j also topples, distributing 
a fraction of stress aFj to its four nearest neighbors. Such a sequence of topplings 
continues until the stress of all sites on the lattice becomes smaller than the threshold 
Fth- A sequence of toppling events, which is assumed to occur instantaneously, cor- 
responds to one seismic event or an avalanche. After an avalanche, the system goes 
into an interseismic period where uniform loading of F is resumed, until some of the 
sites reach the threshold and the next avalanche starts. The transmission parameter a 



2 The Model 




(i) 
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measures the extent of nonconservation of the model. The system is conservative for 
a = 0.25, and is nonconservative for a < 0.25. 

Improvement on original OFC model. It has been widely accepted that earthquakes 
occur as a result of the relative motion of tectonic plates. The plates move relatively 
to one another, resulting in the build up of stress at the plate boundaries. When the 
stress at the plate boundaries reaches to a level that cannot be supported by friction 
between the plates, the strain energy is released intermittently, that is, an earthquake 
happens. We notice that the tectonic plates which bear higher stress are prone to be 
affected by other plates. It can collect more energy or force released by other plates. 
So it is reasonable that the plate with higher stress will get more energy or force when 
its adjacent plate is released. Based on the argument above, we can assume the edge 
weight Wij{t) — [Fi(t) + Fj(t)] /2, for the simplicity of our model, which is deter- 
mined by the seismogenic forces of the two connected sites. This assumption is not 
only a good simulation of the above, but also more importantly, it can be used to model 
the heterogeneity of energy transfer. In order to study the dynamics of our weighted 
OFC model, we should reconsider the redistribution rule. Compared with the original 
OFC model, we just need a new transmission parameter ctj defined as below 11311201 : 



otjit) = a x ^ — = a x — ^ , „ — — ' , „ . (2) 

E 4 *;•(«)+ £ F^t) 



In our improved model, the factor ctj (t) defines the level of local conservation of the 
system and can be adjusted by parameter a. Therefore, for the sake of convenience, 
we consider the parameter a as the control parameter. For a generic initial condition, 
the weighted OFC model, after some transients (discarded in each run), builds up long- 
range spatial correlations, reaches a critical state and generate a time series of avalanche 
size Si, i = 1, n. In particular, we will analyze a time series of n — 10 7 events. 

3 Simulation Results 

3.1 Avalanche- Size Distributions-Effect of the control parameter 

In this section, we mainly analyze the probability distribution of the avalanche sizes. 
The weighted OFC model generates an avalanche size sequence and the avalanche size 
distribution is the frequency of the occurrence of the avalanches with the same size. In 
our model, there are a number of adjustable parameters. For example, we can adjust 
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the threshold for each node, so nodes can be considered as special (this issue will be 
given in a future work) or one can also consider the impact of network structure (this 
issue will be addressed below). In this section, we consider the relationship between 
avalanche size distribution and the control parameter a. As in original OFC model, 
the control parameter can also be used to measure avalanche behavior. The continuous, 
nonconservative weighted OFC model exhibits SOC behavior for a wide range of a 
values. The avalanche size exponent r depends on a |20l . 




io 4 a75 

Figure 1 : Distribution of avalanche sizes for different control parameter a, L = 64 and 

Padd = 0. 



In Fig. 1 we plot the avalanche size distribution for the weighted OFC model with 
different control parameters a 1211 . From this figure, we find that the system develops 
an approximate power-law distribution for avalanche sizes in a wide range of parameter 
a. When a < 0.88, the model only produces a power-law distribution for avalanche 
sizes. The system not only shows power-law behavior but also satisfies the finite-size 
scaling in the parameter range a — 0.88 to a = 1. 
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3.2 Avalanche-Size Distribution-finite-size scaling 

To verify the criticality of our weighted model, we study the effect of increasing the 
system size L. We observe that, for each constant value of a, the avalanche size ex- 
ponent t does not change, while the cutoff in the energy distribution scales with the 
system size. Our weighted OFC model not only shows power-law behavior in the 
avalanche size distribution, but also satisfies the finite-size scaling behavior in the pa- 
rameter range mentioned above. In this part, we propose a simple finite-size scaling 
analysis for the avalanche size distribution of the form 

P(S,L) oc L-Pg(S/L u ) , (3) 

where g is the so-called universal scaling function, parameters f3 and v are critical ex- 
ponents used to characterize scaling properties, v may reflect the scaling relationship 
between the cut off of the distribution function and the system size, while (3 is a normal- 
ization parameter. Fig. 2 displays P(S, L) versus the avalanche size S for the weighted 
OFC model on square lattice of size L = 32, 48, 64 with control parameter a = 1 and 
the inset of Fig. 2 displays the transformed avalanche size distribution, L l3 P(S, L), 
versus rescaled avalanche size, S/ L v . A clear data collapse is evident for the proposed 
scaling function with /3 = 2.456, v = 2.002. The value of critical avalanche size ex- 
ponent (r = 1.220 ± 0.003) [20 1 is in agreement with the finite-size scaling hypothesis 
since for asymptotically large N, it is well-known that P(S) ~ S~ T with t — f3/v 
ll22l . which gives r = 1.227 for the obtained values of j3 and v. So far, we can conclude 
that our weighted OFC model are not only self-organized but also critical. 

3.3 Avalanche- Size Distribution-Effect of long range parameter 

In this section, we mainly discuss the effects of long-range connectivity. The reasons 
why we introduce long range to our model must be given first. To begin with, con- 
structing networks from real seismic data, Baiesi and Paczuski as well as Abe and 
Suzuki reported the discoveries of the scale-free and the small-world features in real 
earthquakes ||8] [17J |T8l. Then, according to the geophysics and geology, heteroge- 
neous character of real earthquake systems and the effect of long-range interactions in 
real earthquakes were found by Mori and Kawamura [23 1, for instance in earthquake 
triggering and interaction, where the static stress may involve relaxation processes in 
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Figure 2: The avalanche size distribution P(S, L) for the weighted OFC model on 
square lattice with system size L = 32, 48, 64. In the inset, the transformed avalanche 
size distribution versus rescaled avalanche size is given. 

the asthenosphere with relevant spatial and temporal long-range effects. Here, we in- 
troduce a small fraction of long-range links (denoted as the long range parameter p a dd) 
in the lattice so as to obtain a small world topology. The long-range connections largely 
reduce the average distance of the original network (here, our model is based on the 
NW small- world model). 

In ll20l . the effects of the control parameter have mainly discussed, here we only 
discuss the behavior of the critical state. In this model, the state of the system is con- 
trolled by the control parameter and long rang parameter p a dd- So, whether the system 
is in self-organized criticality or not depends on these two parameters. Depending on 
the long-range parameter, the network can produce a rich repertoire of behaviors. In 
Fig. 3, we fix a — 1 and show examples of avalanche size distributions for various val- 
ues of Padd- For small values of p a dd (Padd < 0.3), critical avalanche size distributions 
are observed. This regime is characterized by an approximate power-law distribution 



7 



for avalanche sizes almost up to the system sizes where an exponential cutoff is ob- 
served. For larger values of p add (p a dd > 0.3), the distribution is supercritical, that is, 
a substantial fraction of triggering events spread through the whole system [21 1. When 
the control parameter equals to other values, the system shows self-organized criticality 
behavior which is different from the behavior above, and this can be explained as it is 
on the more susceptive critical state (a = 1) than others. In the inset of Fig. 3, we have 
simulated this behavior based on different lattices and found that they show the same 
behavior independent of lattice size expect for L < 16, which can be considered as the 
effects of boundary conditions. 
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Figure 3: Avalanche size distributions for different long range parameter p a dd- The 
long range parameter near p a dd < 0.3 seems critical. For p a dd > 0.3, the distribution 
is supercritical. The inset shows P(S, L) versus the avalanche size S for the weighted 
OFC model on square lattice with system sizes L = 10, 16, 32 and 48 for a = 0.97 and 

Padd = 0.7. 
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3.4 Probability density function for the avalanche size difference 

In recent years, SOC models have been intensively studied considering time intervals 
between avalanches in the critical regime ||24| . Here, we follow a different approach 
which reveals interesting information on the eventual criticality of the model under 
examination. Inspired by recent studies on turbulence and the time-series of real earth- 
quakes, we introduce the distribution of returns, i.e., the differences between fluctuation 
lengths obtained at consecutive time steps, as AS(t) — S(t + S) — S(t), on the differ- 
ences between avalanche sizes calculated at time t + 5 and at time t, 5 being a discrete 
time interval ]25l |26l [9) . It should also be noted that, in order to have zero mean, the 
returns are normalized by introducing the variable x as: 



where < . > stays for the mean value of the given data set II27I . The signal of the 
distribution of returns reveals very interesting results on the criticality of the weighted 
OFC model. In recent works [9 , 28 , 29 ], it is shown that the return distributions can be 
well approximated by a g-Gaussian of type 



(which are the standard distributions obtained in nonextensive statistical mechanics 
ll30l [3T | and from where the standard Gaussian form is obtained as a special case for 
q — > 1) when the avalanche size distribution is a power-law with an exponent r. More- 
over, it is also found that the appropriate q value could be determined a priori from the 
exact relation 



given in (29). It is clear from those efforts that, as the system size increases, the power- 
law regime in avalanche size distribution persists more and more (before arriving the 
exponential decay part) which makes the appropriate (/-Gaussian for the return distri- 
bution to dominate more and more the tails together with the central part. On the other 
hand, usually it is very difficult (if not impossible) to reach very large system sizes (in 
order to approach thermodynamic limit) in such model systems. For the small system 
sizes, considerably short power-law regime is immediately followed by the exponential 



x = AS- <S> 



(4) 




(5) 
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decay in avalanche size distribution and consequently this yields in the return distribu- 
tion the appropriate q-Gaussian to deteriorate in the tails l|29ll28ll . 

In order to explain this tendency, a mathematical simple model for finite-size effects 
exhibiting the gradual approach to g-Gaussians, has been proposed using the following 
differential equation PTl[32ll : 

d ' ! . hry r _ { p q _ hr)y i (b q >b r >0;q>l;y(0) = l). (7) 



d(x 2 ) 

For the particular case r = 1, if one takes bi = 0, the solution is given by the q- 
Gaussian y = [l — (1 — q)b q x 2 ] q ^ = e q bqX . If b q — b%, the solution is 
given by the Gaussian y = e~ hl x . For the case b q > bi > and q > 1, we obtain a 
crossover between these two solutions, the |x| — > oo asymptotic one being the Gaussian 
behavior. For this particular case with q > 1, the solution can be found as an explicit 
expression of the form y(x), namely, 

y = - — ■ (8) 

1 _ h. _|_ h. e (g-i)6i x 2 

bi ~ bi 

On the other hand, for the particular case r = with q > 1, the solution can only 
be given by the explicit x(y) form, namely, 

^_»uri,,, 1+ i,_(^i_^ri, 1 , 1+ i,_&zWj,}, (B 

b I* Lg q b i lq q b a J J 

where 2F1 is the hypergeometric function. 

Indeed, the weighted OFC model studied here constitutes a very good example to 
check the validity of both solutions since (i) it is an example of a model which can only 
be simulated with small system sizes and (ii) the model allows us to define two types 
of avalanche definition, one of which seems to produce return distributions that can 
be approached by Eq.®, whereas the other definition yields return distributions that 
can be given by Eq.©. The standard way of defining an avalanche (the one also used 
throughout this work) is to include each triggered site only once during an avalanche 
which restricts the size of an avalanche with the size of the system. This definition 
results in the return distributions shown in Fig. [4] It is easily seen that the returns 
are well approximated by the crossover formula given in Eq.© with q = 2.64 which 
comes a priori from Eq.©. It is also evident from this figure that, as L — > 00 limit 
is approached, it seems that the return distributions would converge to the g-Gaussian 
with q = 2.64 for the entire region including the central part and the tails. 
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Figure 4: The probability distribution functions of the weighted OFC model with re- 
stricted avalanche definition for representative system sizes. The crossover formula 
given in Eq.® seems to describe the tendency in the entire region except the turn- 
ing points in the tails (central part is given in the Inset). As system size increases, it 
is clearly seen that the return distributions appear to approach the prefect ^-Gaussian 
curve better and better. 

Another way of defining an avalanche is to relax the restriction that allows each site 
to trigger only once during an avalanche. This means that, during a running avalanche, 
one site can be triggered more than once which clearly relaxes the restriction of hav- 
ing maximum avalanche sizes of the order of system size. The use of such definition 
does not change the value of the avalanche size distribution exponent r but results 
in a smoother crossover from the power-law regime to exponential decay part. This 
observed tendency would be expected to have an effect also in the return distributions. 
This can be seen in Fig. |5]where the return distributions can now be well approached by 
the crossover formula given in Eq.®. The gradual approach to the perfect g-Gaussian 
is evident as the system size is increased. It is also worth noting here that the ob- 
served behavior of return distributions do not depend on the interval S considered for 
the avalanche size difference, which is also the case for the worldwide and the northern 
California catalogs (9j. 
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Figure 5: The probability distribution functions of the weighted OFC model with non- 
restricted avalanche definition for representative system sizes. The crossover formula 
given in Eq.® seems to perfectly describe the tendency in the entire region (central 
part is given in the Inset). As system size increases, it is clearly seen that the return 
distributions appear to approach the perfect gr-Gaussian curve better and better. 

As a result, one can conclude here that the behavior of the return distributions is 
very different from a Gaussian shape and seems to be well approached by one of the 
two crossover formulas (either by Eq.© if the non-restricted definition of avalanche is 
used or by Eq.© if the restricted definition of avalanche is used). As far as we know, 
this constitutes the first example in literature where the two forms of these crossover 
solutions can be used together in the same model system. Finally, all the numerical 
findings obtained here suggest that, as the thermodynamic limit is approached, the 
behavior of the return distributions seems to converge to the appropriate ^-Gaussian 
shape in the entire region. 

3.5 Statistics of waiting time in weighted OFC model 

Recently, a new necessary SOC signature has been proposed, in the context of solar 
flare dynamics. It is based on a different type of statistics that deals with waiting times, 
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i.e., the time intervals between two successive bursts or avalanches. In this section, we 
consider the waiting time of avalanche sequences generated by our weighted model. 
The waiting time is defined as the time between the first trigger and the second one. We 
use the overall statistical method and statistics of the waiting times for any node, then 
the distribution of all nodes. After this, we calculated probability distribution of waiting 
times. It can be argued that, if the triggers are not correlated, the process should be 
somehow related to a Poisson process, and the probability distribution function of the 
waiting times should be an exponential law. However,the existence of extended power 
laws in the waiting-time probability distribution function of solar flare measurements 
has been noticed by several authors ll33~l [341 . Several years ago, Christensen et al. 
showed that waiting times would follow power-law distributions if only events larger 
than a certain size are considered in the context of a spring -block model for earthquakes 
ll35l . In this paper, we also analyze the waiting time distribution of the weighted OFC 
model which exhibits a clear nonexponential behavior as can be seen in Fig. 6. The 
power-law regime of the waiting time distribution lasts about two decades. 
We propose a scaling relation for the waiting-time distribution of the form 

P(T) oc L~ e g{T / LP) , (10) 

with the scaling exponents 8 = 3.80 and 7 = 1.75, which are shown (see the inset of 
Fig. 6) to be consistent with the data coming from our model. 

4 Summary and Conclusion 

In order to obtain the inhomogeneous network and different local friction and elasticity, 
we have introduced the weighted edge to improve the original redistribution rule. We 
have shown self-organized criticality in the weighted coupled map lattice. The proba- 
bility density functions of the avalanche size differences (namely, return distributions) 
appear to exhibit fat tails that can be approached by a g-Gaussian shape, in the thermo- 
dynamic limit, with an appropriate value of q coming a priori from the avalanche size 
exponent r. Moreover, for the small system sizes, the observed behavior of the return 
distributions seems to obey the crossover formula proposed in ll32l in order to explain 
the transition from the q-Gaussian behavior to the Gaussian observed so far in some 
other model systems with small system sizes. These results could be interpreted that 
there are no correlations between any two seismic behavior. Our findings support the 
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Figure 6: The waiting time distribution and its scaling analysis. 

hypothesis that even the statistical data of previous earthquake is known, the magni- 
tude of the next earthquake is still unpredictable. Finally, the scaling relation of waiting 
times for the weighted OFC model has been discussed and obtained. 
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